#include "VectorOps.h"
#include <iostream>
#include <stdexcept>


using namespace std;


VectorOps::VectorOps(vector<SeamStress::Seamstress*>* ss) : pins(this, ss)
{
  nthreads = ss->size();
}


VectorOps::~VectorOps()
{
  
}


void VectorOps::copyVector(const refloat_vec* input, refloat_vec* output, unsigned long int offset, unsigned long int offset2)
{
  thread_var1 = offset;
  thread_var2 = offset2;
  vecpointer1 = input;
  vecpointer2 = output;
  pins.sewStraight(&VectorOps::copyVectorThread, nthreads);
}


void VectorOps::copyVectorThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  unsigned long int size = vecpointer1->size();
  if( (size - 1 + thread_var2 - thread_var1) >= (vecpointer2->size())  ){size = (vecpointer2->size() - thread_var2 + thread_var1);}
  if( (vecpointer2->size() < vecpointer1->size()) && (thread_var1==0) && (thread_var2==0) ){size = vecpointer2->size();}
  for(unsigned long int i=(w+thread_var1);i<size;i+=nthreads)
  {
    (*vecpointer2)[i + thread_var2 - thread_var1].set((*vecpointer1)[i]);
  }
}


void VectorOps::copyVector(const cxfloat_vec* input, cxfloat_vec* output, long unsigned int offset, long unsigned int offset2)
{
  thread_var1 = offset;
  thread_var2 = offset2;
  cxvecpointer1 = input;
  cxvecpointer2 = output;
  pins.sewStraight(&VectorOps::copyVectorCxThread, nthreads);
}


void VectorOps::copyVectorCxThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  unsigned long int size = cxvecpointer1->size();
  if( (size - 1 + thread_var2 - thread_var1) >= (cxvecpointer2->size())  ){size = (cxvecpointer2->size() - thread_var2 + thread_var1);}
  if( (cxvecpointer2->size() < cxvecpointer1->size()) && (thread_var1==0) && (thread_var2==0) ){size = cxvecpointer2->size();}
  for(unsigned long int i=(w+thread_var1);i<size;i+=nthreads)
  {
    (*cxvecpointer2)[i + thread_var2 - thread_var1].set((*cxvecpointer1)[i]);
  }
}


void VectorOps::copyVector(const refloat_vec* input, cxfloat_vec* output, long unsigned int offset, long unsigned int offset2)
{
  thread_var1 = offset;
  thread_var2 = offset2;
  vecpointer1 = input;
  cxvecpointer2 = output;
  pins.sewStraight(&VectorOps::copyVectorReCxThread, nthreads);
}


void VectorOps::copyVectorReCxThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  unsigned long int size = vecpointer1->size();
  if( (size - 1 + thread_var2 - thread_var1) >= (cxvecpointer2->size())  ){size = (cxvecpointer2->size() - thread_var2 + thread_var1);}
  if( (cxvecpointer2->size() < vecpointer1->size()) && (thread_var1==0) && (thread_var2==0) ){size = cxvecpointer2->size();}
  refloat zero(0, mpfr_get_prec(refloat_cache1->value) );
  for(unsigned long int i=(w+thread_var1);i<size;i+=nthreads)
  {
    (*cxvecpointer2)[i + thread_var2 - thread_var1].set((*vecpointer1)[i], zero);
  }
}


void VectorOps::reverseCopy(const refloat_vec* input, refloat_vec* output)
{
  vecpointer1 = input;
  vecpointer2 = output;
  pins.sewStraight(&VectorOps::reverseCopyThread, nthreads);
}


void VectorOps::reverseCopyThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<vecpointer1->size();i+=nthreads)
  {
    (*vecpointer2)[i].set((*vecpointer1)[vecpointer1->size() - 1 - i]);
  }
}


void VectorOps::reverseCopy(const cxfloat_vec* input, cxfloat_vec* output)
{
  cxvecpointer1 = input;
  cxvecpointer2 = output;
  pins.sewStraight(&VectorOps::reverseCopyCxThread, nthreads);
}


void VectorOps::reverseCopyCxThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<cxvecpointer1->size();i+=nthreads)
  {
    (*cxvecpointer2)[i].set((*cxvecpointer1)[cxvecpointer1->size() - 1 - i]);
  }
}


void VectorOps::negateVector(refloat_vec* input, refloat_vec* output)
{
  vecpointer1 = input;
  vecpointer2 = output;
  pins.sewStraight(&VectorOps::negateVectorThread, nthreads);
}


void VectorOps::negateVectorThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<vecpointer1->size();i+=nthreads)
  {
    (*vecpointer2)[i].set((*vecpointer1)[i]);
    (*vecpointer2)[i].neg();
  }
}


void VectorOps::negateVector(cxfloat_vec* input, cxfloat_vec* output)
{
  cxvecpointer1 = input;
  cxvecpointer2 = output;
  pins.sewStraight(&VectorOps::negateVectorCxThread, nthreads);
}


void VectorOps::negateVectorCxThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<cxvecpointer1->size();i+=nthreads)
  {
    (*cxvecpointer2)[i].set((*cxvecpointer1)[i]);
    (*cxvecpointer2)[i].neg();
  }
}


void VectorOps::conjugateVector(const cxfloat_vec* input, cxfloat_vec* output)
{
  cxvecpointer1 = input;
  cxvecpointer2 = output;
  pins.sewStraight(&VectorOps::conjugateVectorCxThread, nthreads);
}


void VectorOps::conjugateVectorCxThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<cxvecpointer1->size();i+=nthreads)
  {
    const cxfloat& val1 = (*cxvecpointer1)[i];
    cxfloat& val2 = (*cxvecpointer2)[i];
    
    val2.set(val1);
    val2.conj();
  }
}


void VectorOps::realVector(const cxfloat_vec* input, refloat_vec* output)
{
  cxvecpointer1 = input;
  vecpointer2 = output;
  pins.sewStraight(&VectorOps::realVectorThread, nthreads);
}


void VectorOps::realVectorThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<cxvecpointer1->size();i+=nthreads)
  {
    const cxfloat& val1 = (*cxvecpointer1)[i];
    refloat& val2 = (*vecpointer2)[i];
    
    val1.get_real(val2);
  }
}

void VectorOps::reciprocateVector(refloat_vec* input, refloat_vec* output)
{
  vecpointer1 = input;
  vecpointer2 = output;
  pins.sewStraight(&VectorOps::reciprocateVectorThread, nthreads);
}


void VectorOps::reciprocateVectorThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<vecpointer1->size();i+=nthreads)
  {
    (*vecpointer2)[i].set((*vecpointer1)[i]);
    (*vecpointer2)[i].reciprocal();
  }
}


void VectorOps::reciprocateVector(cxfloat_vec* input, cxfloat_vec* output)
{
  cxvecpointer1 = input;
  cxvecpointer2 = output;
  pins.sewStraight(&VectorOps::reciprocateVectorCxThread, nthreads);
}


void VectorOps::reciprocateVectorCxThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<cxvecpointer1->size();i+=nthreads)
  {
    (*cxvecpointer2)[i].set((*cxvecpointer1)[i]);
    (*cxvecpointer2)[i].reciprocal();
  }
}


//output += input
void VectorOps::addVector(const refloat_vec* input, refloat_vec* output)
{
  vecpointer1 = input;
  vecpointer2 = output;
  pins.sewStraight(&VectorOps::addVectorThread, nthreads);
}


void VectorOps::addVectorThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  unsigned long int size = vecpointer1->size();
  if( (vecpointer2->size() < vecpointer1->size()) ){size = vecpointer2->size();}
  
  for(unsigned long int i=w;i<size;i+=nthreads)
  {
    (*vecpointer2)[i].add((*vecpointer1)[i]);
  }
}


void VectorOps::addVector(cxfloat_vec* input, cxfloat_vec* output)
{
  cxvecpointer1 = input;
  cxvecpointer2 = output;
  pins.sewStraight(&VectorOps::addVectorCxThread, nthreads);
}


void VectorOps::addVectorCxThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  unsigned long int size = cxvecpointer1->size();
  if( (cxvecpointer2->size() < cxvecpointer1->size()) ){size = cxvecpointer2->size();}
  
  for(unsigned long int i=w;i<size;i+=nthreads)
  {
    (*cxvecpointer2)[i].add((*cxvecpointer1)[i]);
  }
}


void VectorOps::subVector(cxfloat_vec* input, cxfloat_vec* output)
{
  cxvecpointer1 = input;
  cxvecpointer2 = output;
  pins.sewStraight(&VectorOps::subVectorCxThread, nthreads);
}


void VectorOps::subVectorCxThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  unsigned long int size = cxvecpointer1->size();
  if( (cxvecpointer2->size() < cxvecpointer1->size()) ){size = cxvecpointer2->size();}
  
  for(unsigned long int i=w;i<size;i+=nthreads)
  {
    (*cxvecpointer2)[i].sub((*cxvecpointer1)[i]);
  }
}


// output = input1 + scale*input2
void VectorOps::addMul(const refloat_vec* input1, const refloat_vec* input2, refloat_vec* output, refloat& scale)
{
  if(input1->size() != output->size()){throw std::runtime_error("VectorOps::addMul : input and output must be of the same size");}
  if(input2->size() != output->size()){throw std::runtime_error("VectorOps::addMul : input and output must be of the same size");}
  
  vecpointer1 = input1;
  vecpointer4 = input2;
  vecpointer2 = output;
  val1 = &scale;
  pins.sewStraight(&VectorOps::addMulThread, nthreads);
}


void VectorOps::addMulThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  refloat temp1(*val1);
  
  for(unsigned long int i=w;i<vecpointer1->size();i+=nthreads)
  {
    temp1.set( (*vecpointer4)[i] );temp1.mul(*val1);temp1.add( (*vecpointer1)[i] );
    (*vecpointer2)[i].set(temp1);
  }
}


//output -= input
void VectorOps::subVector(const refloat_vec* input, refloat_vec* output)
{
  try
  {
    if(input->size() != output->size()){throw (char*)("subVector arguments must be of same size");}
  }
  catch(char* str)
  {
    cout<<"Exception from VectorOps: "<<str<<endl;
    throw;
    return;
  }
  vecpointer1 = input;
  vecpointer2 = output;
  pins.sewStraight(&VectorOps::subVectorThread, nthreads);
}


void VectorOps::subVectorThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<vecpointer1->size();i+=nthreads)
  {
    (*vecpointer2)[i].sub((*vecpointer1)[i]);
  }
}


void VectorOps::elementWiseMultiply(const refloat_vec* input, refloat_vec* output)
{
  vecpointer1 = input;
  vecpointer2 = output;
  pins.sewStraight(&VectorOps::elementWiseMultiplyThread, nthreads);
}


void VectorOps::elementWiseMultiplyThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  unsigned long int size = vecpointer1->size();
  if(vecpointer2->size() < vecpointer1->size()){size = vecpointer2->size();}
  
  for(unsigned long int i=w;i<size;i+=nthreads)
  {
    (*vecpointer2)[i].mul((*vecpointer1)[i]);
  }
}


void VectorOps::elementWiseMultiply(const cxfloat_vec* input, cxfloat_vec* output)
{
  cxvecpointer1 = input;
  cxvecpointer2 = output;
  pins.sewStraight(&VectorOps::elementWiseMultiplyCxThread, nthreads);
}


void VectorOps::elementWiseMultiplyCxThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  unsigned long int size = cxvecpointer1->size();
  if(cxvecpointer2->size() < cxvecpointer1->size()){size = cxvecpointer2->size();}
  
  for(unsigned long int i=w;i<size;i+=nthreads)
  {
    (*cxvecpointer2)[i].mul((*cxvecpointer1)[i]);
  }
}


//output /= input
void VectorOps::divVector(refloat_vec* input, refloat_vec* output)
{
  try
  {
    if(input->size() != output->size()){throw (char*)("divVector arguments must be of same size");}
  }
  catch(char* str)
  {
    cout<<"Exception from VectorOps: "<<str<<endl;
    throw;
    return;
  }
  vecpointer1 = input;
  vecpointer2 = output;
  pins.sewStraight(&VectorOps::divVectorThread, nthreads);
}


void VectorOps::divVectorThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<vecpointer1->size();i+=nthreads)
  {
    (*vecpointer2)[i].div((*vecpointer1)[i]);
  }
}


void VectorOps::divVector(cxfloat_vec* input, cxfloat_vec* output)
{
  try
  {
    if(input->size() != output->size()){throw (char*)("divVector arguments must be of same size");}
  }
  catch(char* str)
  {
    cout<<"Exception from VectorOps: "<<str<<endl;
    throw;
    return;
  }
  cxvecpointer1 = input;
  cxvecpointer2 = output;
  pins.sewStraight(&VectorOps::divVectorCxThread, nthreads);
}


void VectorOps::divVectorCxThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<cxvecpointer1->size();i+=nthreads)
  {
    (*cxvecpointer2)[i].div((*cxvecpointer1)[i]);
  }
}


void VectorOps::scaleVector(refloat_vec* input, refloat& scale)
{
  vecpointer2 = input;
  val1 = &scale;
  pins.sewStraight(&VectorOps::scaleVectorThread, nthreads);
}


void VectorOps::scaleVectorThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<vecpointer2->size();i+=nthreads)
  {
    (*vecpointer2)[i].mul(*val1);
  }
}


void VectorOps::scaleVector(cxfloat_vec* input, cxfloat& scale)
{
  cxvecpointer2 = input;
  cxval1 = &scale;
  pins.sewStraight(&VectorOps::scaleVectorCxThread, nthreads);
}


void VectorOps::scaleVectorCxThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<cxvecpointer2->size();i+=nthreads)
  {
    (*cxvecpointer2)[i].mul(*cxval1);
  }
}


void VectorOps::dotProduct(const cxfloat_vec* input1, const cxfloat_vec* input2, refloat& value)
{
  try
  {
    if(input1->size() != input2->size()){throw (char*)("dotProduct inputs must be of same size");}
  }
  catch(char* str)
  {
    cout<<"Exception from VectorOps: "<<str<<endl;
    throw;
    return;
  }
  
  cxvecpointer1 = input1;
  cxvecpointer4 = input2;
  val1 = &value;
  thread_tot = 1;
  unsigned long int tt=1;
  while(tt <= nthreads && tt<=input1->size())
  {
    thread_tot = tt;
    tt <<= 1;
  }
  refloat_vec thread_vec;
  refloat temp(value);temp.set(0);
  thread_vec.assign(thread_tot, temp);
  vecpointer3 = &thread_vec;
  pins.sewStraight(&VectorOps::dotProductCxThread1, thread_tot);
  while(thread_tot >= 2)
  {
    thread_tot >>= 1;
    pins.sewStraight(&VectorOps::dotProductThread2, thread_tot);
  }
  value.set(thread_vec[0]);
}


void VectorOps::dotProductCxThread1(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  refloat temp((*vecpointer3)[w]);
  
  cxfloat temp1((*cxvecpointer4)[0]);
  cxfloat temp2((*cxvecpointer4)[0]);
  
  for(unsigned long int i=w;i<cxvecpointer1->size();i+=thread_tot)
  {
    temp1.set((*cxvecpointer4)[i]);temp1.conj();
    mul((*cxvecpointer1)[i], temp1, temp2);
    temp2.get_real(temp);
    (*vecpointer3)[w].add(temp);
  }
}


void VectorOps::dotProduct(const refloat_vec* input1, const refloat_vec* input2, refloat& value)
{
  try
  {
    if(input1->size() != input2->size()){throw (char*)("dotProduct inputs must be of same size");}
  }
  catch(char* str)
  {
    cout<<"Exception from VectorOps: "<<str<<endl;
    throw;
    return;
  }
  
  vecpointer1 = input1;
  vecpointer4 = input2;
  val1 = &value;
  thread_tot = 1;
  unsigned long int tt=1;
  while(tt <= nthreads && tt<=input1->size())
  {
    thread_tot = tt;
    tt <<= 1;
  }
  refloat_vec thread_vec;
  refloat temp(value);temp.set(0);
  thread_vec.assign(thread_tot, temp);
  vecpointer3 = &thread_vec;
  pins.sewStraight(&VectorOps::dotProductThread1, thread_tot);
  while(thread_tot >= 2)
  {
    thread_tot >>= 1;
    pins.sewStraight(&VectorOps::dotProductThread2, thread_tot);
  }
  value.set(thread_vec[0]);
}


void VectorOps::dotProductThread1(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  refloat temp((*vecpointer3)[w]);
  
  for(unsigned long int i=w;i<vecpointer1->size();i+=thread_tot)
  {
    mul((*vecpointer1)[i], (*vecpointer4)[i], temp);
    (*vecpointer3)[w].add(temp);
  }
}


void VectorOps::dotProductThread2(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  (*vecpointer3)[w].add((*vecpointer3)[w + thread_tot]);
}


void VectorOps::setVector(refloat_vec* input, const refloat& val)
{
  vecpointer2 = input;
  val1 = &val;
  pins.sewStraight(&VectorOps::setVectorThread, nthreads);
}


void VectorOps::setVectorThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<vecpointer2->size();i+=nthreads)
  {
    (*vecpointer2)[i].set(*val1);
  }
}


void VectorOps::setVector(cxfloat_vec* input, const cxfloat& val)
{
  cxvecpointer2 = input;
  cxval1 = &val;
  pins.sewStraight(&VectorOps::setVectorCxThread, nthreads);
}


void VectorOps::setVectorCxThread(void* arg)
{
  unsigned long int w = (*((unsigned long int*)arg));
  
  for(unsigned long int i=w;i<cxvecpointer2->size();i+=nthreads)
  {
    (*cxvecpointer2)[i].set(*cxval1);
  }
}












